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5 Abstract 

We present a new method in order to get variable slip coefficient in binary 
lattice Boltzmann models to simulate gaseous flows. Boundary layer theory is 
presented. We study both the single- and multi-fluid BGK-typc models as well. 
The boundary slip and the Knudscn layer are analyzed in detail. Benchmark 
simulations are carried out in order to compare the analytical derivation with 
the numerical results. Excellent agreement is found between the two situations. 
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8 Introduction 

Over the last decades, the lattice Boltzmann model (LB) has become a versatile 
tool in modeling complex hydrodynamic problems pQ . Recently, it has attracted 
considerable attention in modeling rarefied gas flows. These flows arc important 
in several applications ranging from micro-electromechanical systems (MEMS) 
[5] to aeronautical applications. With increasing rarefaction characterized by the 
Knudsen number, the ratio of the molecular mean- free path and the macroscopic 
sizes, the flow experiences some interesting phenomena, such as the gaseous 
slip along the solid wall. The slip phenomenon is one of the cornerstone of 
micro-channel applications. In LB, several methods have been proposed for 
the description of the slip flow, mostly based on kinetic boundary conditions 
[31 HI [SI IS1 El [HI [SI [IH] , and analytical formulas have been derived for the slip 
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velocity and the slip coefficient. We mention that the first kinetic boundary 
condition was introduced by Lim et al. |10j . The above approaches focused on 
one-component gases and the slip phenomenon in mixture LB models is yet not 
understood. In addition, it is worth mentioning that the bare diffuse reflection 
boundary condition can not reproduce the correct value of the slip coefficient 

In this paper, we present a boundary layer theory for binary lattice Boltz- 
mann models and introduce a new boundary treatment in order to get tunable 
slip coefficient in the method. We consider both single and multicomponent 
BGK-typc collision operators as well. In the framework of the boundary layer 
theory, analytical formula is derived for the slip coefficient. Simulations are per- 
formed to validate the model. The analytical results are in excellent agreement 
with the ones obtained from the numerical simulations. 

9 The binary mixture model 

We begin our consideration with the binary discrete Boltzmann equation for 
single- and multi-fluid BGK-typc collision operators 

d t f?(x,t)+c aia d a f[( X ,t) = K+^ m )[/r (e) (p CT ,M CT )-/n + 

^[r (e) (p CT ,«)-/r (e) (pa,^)], (1) 

where /f,/^ denote the single-specie distribution function and equilibrium 
distribution function, c„i a are the discrete speed vectors and v\ , v m are the 
collision frequencies related to the transport coefficients. 

The macroscopic quantities are obtained as the moments of ff as 

n ° = J2fk' n ° u ° = -fk c <yk, (2) 

k k 
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where n ai u a are the particle density and velocity, respectively. Further, 



p a = n a m a , 



(3) 



where p ai m a are the specie mass density and the specie mass, respectively. The 
mixture macroscopic quantities are defined by 



(4) 



where n, p, u denote the mixture particle density, mixture mass density and 
mass-averaged mixture velocity, respectively. The generic equilibrium distribu- 
tion function is chosen as the second order truncation of the Maxwellian 



/f (e) (p<r> M ) = Wi — 



Pa 



2d 



(5) 



We work in two dimensions, and the discrete speed vectors are given by 
c-ai = c a £i, where is the defined as the D2Q9 non-dimensional velocity model 



(0,0) 



e, = < 



V3cos( - i - -), V3sin(-t 



for i = 

— ) ) for i = l. 

2 ') 

7T . 7T 



("v/Gcosf— i ),V6sm(—i )) for i = 5. 

\ 2 4 2 4 / 



The corresponding weights arc defined by 



(6) 



4/9 for % = 
1/9 for i=l. 
1/36 for i = 5 . 



(7) 



Note that this velocity model is different from the original nine-speed model [13] 
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The velocity vectors are rescaled with the multiplier. This is a convenient 
choice for our purposes and our calculations. The definition of the model is given 
by the above equations, Eq. (|6l7p . To complete our definition, we mention that 
the non-dimensional sound speed of the model is unity. The specie sound speed 
is given by c a = c\Jm/m a , where m = p/n is the averaged mass and c is the 
mixture sound speed. 

Before going to the boundary layer theory, some comments are in order. 
The collision operator in the governing equation, Eq. |T]), is equivalent with 
the so-called Hamel model [H] , which describes coupled relaxation towards the 
local specie and mixture equilibriums [HI [15]. With the choice of v\ = 0, the 
collision model results in the single-fluid BGK model. In the two-fluid model, 
the mixture viscosity is given by p = pc 2 / (y\ + v m ) and the diffusivity takes the 
value of D = m 2 c 2 / ' {rn,\ra2V m ) . 



10 Boundary layer theory 

We turn our attention to the half space problem experienced by a gas flow along 
a solid wall. This is the so-called Kramers problem in kinetic theory [16] . The 
solid wall lies in the y direction and is located at x = 0. In the x > region, the 
gas medium exhibits a shear flow. Under these circumstances, we are looking for 
the macroscopic velocity obtained from Eq. |T]) . This problem can be simplified 
by introducing a new reduced distribution function instead of f" 

F? = m CT c ff V3[/ 6 CT -/f], 
F£ = m a c a V3{fZ-fZ}, 

F£ = m CT c CT V3[/ 5 ff -/ 8 CT ]. (8) 
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We also introduce the following weights, u>i, and one-dimensional velocity vec- 
tors, 6j, corresponding to the new distribution function in such a way that 

wi = l/6, w = 2/3, w 2 = l/6, (9) 
ei = -\/3, e = 0, e 2 = y/3. (10) 

With using the reduced distribution function, the macroscopic specie velocity 
is obtained by PaU ay = ^ fc F% . 

The reduced distribution function obeys the following governing equation 
obtained from Eq. (JTJ) in the steady state 

C<7 eid x F° = v^iJ^Fk +VmUi—^J2 F k ~ fa ( n ) 

k P k o' 

This equation is the boundary layer equation, which needs to be solved for 
the macroscopic velocity. The solution of this differential equation is obtained 
in the form F? = F° exp(A;c). As a consequence, we obtain the following 
generalized eigenvalue problem 

c^xf[ (o) = w e k (o) + v m u,A E E fa + -™)^r (0 U2) 

k P k cr' 

This equation is written out in components in the Appendix, Eq. (|2ip . The 
eigenvalues of the problem are obtained by A = [0, 0, Ao, — Ao], where Ao is listed 
in the Appendix, Eq. (|22p . The two zero eigenvalues correspond to a linear flow 
profile, while the other two to a growing and a decaying exponential solution. We 
remind that in the single component D2Q9 models the non-zero eigenvalues are 
absent. However, these two solutions in the mixture case describe the Knudscn 
layer function fading away in order of the mean free path far from the wall. 
The linear flow profile corresponds to the hydrodynamic shear flow in the 
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bulk region. After some straightforward algebra, one can obtain the linear 
solution from Eq. pip as 

F x = a W 4>_^-_^-_ (13) 

where the H superscript denotes the hydrodynamic part of the distribution 
function and the velocity. 

The general solution of the boundary layer equation, Eq. (jllj) . is given by 
the linear combination of the hydrodynamic shear solution, F a ( H \ and a Knud- 
scn layer part corresponding to the nonzero eigenvalues. Here, we mention that 
in our work, the Knudsen layer concept is identical with the boundary layer 
provided by the established model and obtained from the analytical derivation. 
This is because the LB model is a special discrete ordinate method for the solu- 
tion of the Boltzmann equation p~7l[T8]. As it has been recently shown, Knudsen 
layer appears in discrete kinetic models once the reflected populations have dif- 
ferent Ci X normal components |18j . As a result, the boundary layer appearing 
in the model is physically relevant. However, in order to obtain quantitatively 
correct information, one needs to apply optimized boundary condition presented 
below. 

Before determining the slip coefficient in the model, we theoretically show 
that the shear stress is independent of the Knudsen layer part. Indeed, the shear 
stress obtained from the Knudsen layer solution is given by P^ y = m a c a ^ x c 
Si Scr c a e iFi^ exp(Ax). It is quickly realized that P xy is always zero for A 7^ 0, 
because J2a c <r e j-^j is nothing just the sum of the left hand side of all 
equations of Eq. (|21[) . which is identically zero because the right hand side 
becomes the collision invariant of the momenta for the overall mixture. This 
is an important property of the half space problem and remains valid for other 
velocity models as well. As a result, the shear stress at the wall can be obtained 
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from the hydrodynamic solution, Eq. (|13[) . by 

r> \ ' \ " t,<?{H) C 2 du { y H) 

p " = ^ = - ^+. m fa ■ (14) 

In this way, the velocity gradient is always available from the shear stress inde- 
pendently of the unknown Knudsen layer function. 



11 Variable slip coefficient 

In rarefied gas flows, the slip velocity, u s , is considered as the extrapolated gas 
velocity at the wall without the Knudsen layer. It is given in the following form 

"• = OB -p-fa-' (15) 

where a v , (j,,Vo, P are the slip coefficient, the viscosity, the mixture reference 
speed and the pressure, respectively. In the present model, «o = v2c and 
P = pc 2 . In Eq. Q15p. a v is the so-called viscosity based slip coefficient [T5l 
which is a useful definition, since it connects the slip coefficient to a transport 
coefficient, which is always available independently of the molecular details of 
the rarefied gas. In order to get the slip coefficient, we need to determine the 
magnitude of the hydrodynamic solution in the boundary layer. This is fixed 
by the used boundary condition at the wall. 

We introduce a modified diffuse boundary condition in the model in order 
to get variable slip coefficient. The boundary treatment for the reflected popu- 
lations, ei x > 0, are written by 



ft = /.f (e V,0) + (1 - a ) Wl C -^pL^P xy , (16) 

cr pm 



where a is a yet unknown parameter, which will be determined below. For the 
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quantity, F[ , the above boundary condition is written by 



Ff = (l-a)u> 2 V3^^P xy . (17) 
c p 

The particular solution of Eq. (fTTj) is given by the linear combination of the 
decaying exponential solution and the hydrodynamic part. The magnitude of 
hydrodynamic solution and the Knudsen layer part is determined by the bound- 
ary condition, Eq. (jTTJ) , in such a way that 

[ -t \ fcCi Pi o V3ci du y H) 

(1 - a)w 2 V3- — P xy = uj 2 p\u s - uj 2 p\ ■ 5 — 

cr p fx + v m ox 

+kFl {0 \ (18) 

/, \ p;C2 P2 D V3c 2 du y H) 

(1 - fl W 2 V3^- P xy = UJ 2 p 2 U s - UJ 2 p 2 ; 

cr p vi + v m ox 

+fcF 2 2(0) . (19) 



These equations need to be solved for the unknown quantities, u s ,k. It is 
realized that for obtaining the slip velocity, we need the ratio F^ 1 /F 2 2 ^ corre- 



sponding to the decaying eigenvalue — Ao. This can be obtained from Eq. (|21| 
and the result is given in the Appendix, Eq. (1231) . 

After some straightforward calculation, we obtain the slip velocity from Eq. 
(Umi) as 

IIV OUy 

u 8 =aB— (20) 

where the obtained value of the B parameter is given by Eq. in the Ap- 
pendix. It can be seen that we get the desired value of the slip coefficient with 
the choice a = a v / B. 
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Figure 1: Comparison between the measured slip coefficient obtained 
from the LB simulation versus the theoretical derivation at different val- 
ues of the parameter, a. ; A, ▲ correspond to the values of mi/m 2 = 
[1/9.98, 1/9.98, 1/32.804], m/n = [0.5,0.2,0.6], respectively, and represent the 
results obtained from the LB simulations. 4 denotes the analytical results. The 
collision frequencies are chosen by Vi — l/(0.05(iV — 1)), v m = l/(0.03(iV — 1)). 

12 Simulations 

Simulations were performed in order to validate the analytical results obtained 
above. The governing equation Eq. (TTJ was solved using the finite difference 
method. We mention that the finite difference realization is a useful choice 
because of the independence of the velocity vectors and the coordinate grid 
[2T| [22] . The time and the space derivatives are computed using the modified 
Rungc-Kutta and midpoint method, respectively [21]. Note that this scheme 
is second order accurate in both time and space [23j . The Kramers problem 
is simulated in an N x 1 domain, N = 128. The left hand side wall is at 
rest, while the right hand the side moves with a constant velocity to maintain 
the constant shear flow in the domain. The proposed boundary condition with 
the variable slip coefficient is applied at the left wall, where the slip coefficient 
is measured through the relation of Eq. (| 1 5[) . Fig. 1 presents the measured 
slip coefficient versus the theoretical value obtained from Eq. (|2"P)) at different 
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values of a. Excellent agreement is obtained between the two situations. The 
proposed model can be used to tune the value of the slip coefficient. Note that 
above boundary condition can be used in non-stationary case as well. 

Before concluding our work, we mention that the boundary layer theory 
can be further developed involving other flow problems, such as pressure driven 
flows. In that case, the boundary layer equation can be also solved which pro- 
vides a useful background to find an optimal boundary treatment. An another 
interesting issue is the description of complex geometries, slip flow on curved 
surfaces. For the latter situation, Ref. [55] can be generalized, which describes 
variable slip coefficient in the single component LB model on curved boundaries. 

13 Conclusion 

In this paper, we have developed a new boundary treatment in the LB method 
for binary mixtures in order to obtain variable slip coefficient. We have set up a 
boundary layer analysis and derived the analytical value of the slip coefficient. 
Computer simulations have been performed to validate the model. The results 
of the simulations arc in excellent agreement with the analytical derivation. Our 
method can be used to obtain the desirable value of the slip coefficient in the 
lattice Boltzmann model. 
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14 Appendix 



The generalized eigenvalue problem, Eq. (|12p , written out in components is 
given by 



1(0) 



6^ - ™6p 



6^ fc ' "*6 p 

fc r k 



1(0) 
2 l 



2(0) 



6^ • ™6p 



o = < E ^ (0) + <f E(^ 1(0) + *f >) - + ^ 2( 



6 V fc 6 ^ 



2 2(0) - 



c^AF^ = ,4E^ (0) +^4?E(^ (0) +^ (0) )-to+^ 

(21) 



The non-zero eigenvalue of the eigenvalue problem is obtained by 



^/u^ \J m l m 2 V\ + I'm , nnS 

Ao = — , = . \A^> 

■sjv\ + 'iv m m c 

The ratio of the eigenvalues, F^^/F^ ^, corresponding to the decaying 
eigenvalue — Aq can be obtain from Eq. (1211) in such a way that 



F^ = c%(vx + v m + ciV3Xo)k, 

F 2 2(0) = -c\{vi + v m + c 2 V3X )k, (23) 
where k is an arbitrary multiplier. 
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The parameter in the slip velocity, Eq. (|20p , is given in such a way that 
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